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Abstract 


Many  flows  of  aeronautical  interest  have  regions  where  turbulence  has  a  significant  effect.  For 
many  of  these  flows,  Reynolds-averaged  Navier-Stokes  simulation  (RANSS)  techniques  do  not 
give  an  acceptable  description  of  the  flow.  In  these  cases  a  more  detailed  simulation  of  the 
turbulence  is  required.  One  such  detailed  simulation  technique,  large-eddy  simulation  (LES) 
has  matured  to  the  point  of  application  to  complex  flows.  Historically,  LES  have  been  carried 
out  with  structured  grids  which  suffer  from  two  major  difficulties:  the  extension  to  higher 
Reynolds  numbers  leads  to  an  impractical  number  of  grid  points,  and  most  real  world  flows 
are  rather  difficult  to  represent  geometrically  with  structured  grids.  Unstructured-grid  methods 
offer  a  release  from  both  of  these  constraints.  Within  this  sponsored  research  significant  progress 
has  been  made  towards  the  application  of  the  above  approach  to  flows  of  aeronautical  interest. 


1  Introduction 


Computation  of  turbulence  can  be  carried  out  in  three  distinct  ways:  direct  numerical  simula¬ 
tion  (DNS),  where  the  computational  method  resolves  all  of  the  turbulent  motions  (from  the 
largest  scale  down  to  the  scale  where  motion  is  converted  to  heat  via  viscous  dissipation),  LES, 
where  the  large  eddies  carrying  most  of  the  energy  are  resolved  by  the  computational  method 
(leaving  the  small  or  subgrid-scale  (SGS)  eddies  to  be  modeled),  and  RANSS,  where  no  attempt 
is  made  to  resolve  any  of  the  turbulent  motions  (rather,  the  net  effect  of  these  motions  upon 
the  mean  flow  is  modeled).  DNS  can  be  expected  to  remain  out  of  reach  of  most  engineering 
problems  for  several  years.  RANSS,  though  computationally  efficient  for  some  steady  flows,  has 
been  disappointing  for  unsteady  and/or  complex  flows.  Therefore,  LES,  has  been  identified  as 
the  most  promising  turbulence  simulation  tool  for  the  near  future. 

Over  the  past  three  decades,  a  significant  amount  of  research  has  been  devoted  to  LES.  The 
overwhelming  majority  of  this  research  has  been  carried  out  with  spectral  or  structured  grid 
finite  difference  methods.  While  much  has  been  learned  about  the  basic  physics  of  turbulence 
using  these  approaches,  they  suffer  from  two  significant  problems:  difficulty  in  application 
to  arbitrarily  complex  geometries,  and  a  rapid  increase  in  the  number  of  degrees  of  freedom 
necessary  to  describe  flows  at  typical  engineering  Reynolds  numbers.  Recently,  LES  has  been 
extended  to  finite  element  methods  on  unstructured  grids  ([1],[2])  with  near  perfect  parallel 
scaling  on  a  variety  of  different  computing  platforms  [3].  This  extension  not  only  facilitates 
simulation  of  flows  within  or  around  complex  geometries  found  in  engineering  applications,  but 
also  allows  great  reductions  in  computational  effort  through  the  ability  of  unstructured  grids 
to  adapt  locally  to  fine  scale  structures  in  one  flow  region  while  remaining  coarse  in  regions 
.  where  the  structures  are  large.  Already,  these  methods  have  realized  a  27  fold  point  reduction 
over  an  equivalent-resolution  structured  grid  (i.e.  to  make  a  structured  grid  which  matches  the 
resolution  in  the  finest  regions  of  the  unstructured  grid  would  require  a  factor  of  27  times  as 
many  points  assuming  the  stretched  mesh  can  be  mapped  to  a  cube  ([4], [5])). 

2  Method  and  advances  under  AFOSR  support 

2.1  Finite  element  formulation 

There  are  many  finite  element  and  finite  volume  formulations  available  for  application  with 
unstructured  grids.  However,  few  have  had  their  accuracy  and  stability  as  carefully  scrutinized 
as  what  have  become  known  as  stabilized  finite  element  methods.  Two  important  members  of 
this  family,  Galerkin/Least-Squares  (GLS)  and  Streamline  Upwind  Petrov  Galerkin  (SUPG), 
have  been  proven  stable  and  higher  order  accurate  (converging  at  the  optimal  rate  for  a  given 
function  space)  for  flows  ranging  from  inviscid  to  viscous  dominated.  The  early  analyses  were 
done  in  the  context  of  multi-dimensional  linear  advective-diffusive  systems  (see  Hughes  et  al.  [6] 
and  Hughes  et  al.  [7])  and  later  followed  by  an  analysis  of  the  frozen  coefficient  Navier-Stokes 
equations  (see  Franca  and  Hughes  [8]).  The  interested  reader  is  encouraged  to  look  to  these 
references  for  detail  as  the  space  of  this  report  allows  only  a  summary  of  the  formulation. 

Consider  the  compressible  Navier-Stokes  equations  (complete  with  continuity  and  total  en¬ 
ergy  equations)  written  in  filtered  form  after  the  application  of  a  subgrid-scale  model  (see  Moin 
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et  al.  [9],  Moin  &  Jimenez  [10],  Germano  et  al.  [11]  or  section  2.3). 
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Here,  we  use  the  overbar  to  denote  an  unweighted  filter  and  a  tilde  to  denote  a  density  weighted 
filter.  The  filtered  variables  are:  the  velocity  U{,  the  pressure  p,  the  density  p,  the  temperature 
T  and  the  total  energy  et0t-  The  constitutive  laws  relate  the  stress,  r^,  to  the  deviatoric 
portion  of  the  strain,  Sfj  =  «%  —  | Skk$ij,  through  a  molecular  viscosity,  p,  plus  turbulent 
viscosity,  pr-  Similarly,  the  heat  flux,  </;,  is  proportional  to  the  gradient  of  temperature  with  the 
proportionality  constant  given  by  the  addition  of  a  molecular  conductivity,  k,  and  a  turbulent 
conductivity,  nT,  which  is  assumed  proportional  to  the  turbulent  viscosity  as  described  above. 
While  the  formulation  is  not  limited  to  an  ideal  gas,  p  =  pRT,  and  constant  specific  heats  at 
constant  pressure,  cp,  and  at  constant  volume,  c„,  that  is  the  model  used  in  the  calculations 
shown  in  this  report.  Furthermore,  since  all  the  calculations  shown  in  this  report  are  at  low 
Mach  number  where  temperature  variation  is  low  we  have  also  assumed  a  constant  molecular 
viscosity  and  constant  conductivity  through  a  constant  Prandtl  number,  though,  this  again  is 
not  a  necessary  simplification.  Finally  S  is  a  body  force  (or  source)  term. 

For  the  specification  of  the  methods  that  follow,  it  is  helpful  to  define  the  quasi-linear 
operator  (with  respect  to  some  yet  to  be  determined  variable  vector  Y )  related  to  (1)  as 

-  _  .  d  d  3,9.  /r-\ 

c  =  A°m + -  (5) 

from  which  C  can  be  naturally  decomposed  into  time,  advective,  and  diffusive  portions 

C,  =  C,\  +  £a dv  T  (6) 


Here  At  =  F^y  *s  t^ie  ^  Euler  Jacobian  matrix,  Kij  is  the  diffusivity  matrix,  defined  such  that 
KijYj  =  Ffs,  and  A0  =  U  Y  is  the  change  of  variables  metric.  For  a  complete  description  of 
A0,Ai  and  K^,  the  reader  is  referred  to  [12].  Using  this,  we  can  write  (1)  as  simply  CY  =  S. 

To  proceed  with  the  finite  element  discretization  of  the  Navier-Stokes  equations  (1),  we 
must  define  some  notation.  First,  let  12  represent  the  closure  of  the  physical  spatial  domain  (i.e. 
12  U  T  where  T  is  the  boundary).  Next,  12,  is  discretized  into  nei  finite  elements,  12®.  To  derive 
the  weak  form  of  (1),  the  entire  equation  is  dotted  from  the  left  by  a  vector  of  weight  functions, 
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W,  and  integrated  over  the  spatial  domain.  Integration  by  parts  is  then  performed  to  move 
the  spatial  derivatives  onto  the  weight  functions  (reducing  the  continuity  requirements).  This 
process  leads  to  the  integral  equation  (often  referred  to  as  the  weak  form):  find  Y  such  that 

0  =  f  (W  ■  U,t  -  Wti  ■  Ffv  +  Wti  •  Ff ff)  dSl-  j^W-  (-Ffv  +  Ff I)  m  dT 

nel  /» 

+  £  /  CTW  •  r  (CY  -  S )  dfi  (7) 

e=l 

The  first  line  of  (7)  contains  the  Galerkin  approximation  (interior  and  boundary)  and  the  second 
line  contains  the  least-squares  stabilization.  SUPG  stabilization  is  obtained  by  replacing  CT 
by  £adv-  The  stabilization  matrix  r  is  an  important  ingredient  in  these  methods  and  is  well 
documented  in  Shakib  [13]  and  in  Franca  and  Frey  [14].  Note  that  we  have  chosen  to  find  Y 
instead  of  U.  As  discussed  in  Hauke  and  Hughes  [15],  U  is  often  not  the  best  choice  of  solution 
variables,  particularly  when  the  flow  is  nearly  incompressible.  For  the  calculations  performed 
herein,  the  SUPG  stabilized  method  was  applied  with  linearly  interpolated  pressure-primitive 
variables  viz. 

Yi  V 

Y2  ui 
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By  inspecting  (2)-(4)  it  is  clear  that  all  quantities  appearing  in  (7)  may  be  easily  calculated 
from  (8). 

To  develop  a  numerical  method,  the  weight  functions  (W),  the  solution  variable  (Y),  and 
its  time  derivative  (Y>t)  are  expanded  in  terms  of  basis  functions  (typically  piecewise  polyno¬ 
mials;  all  calculations  described  herein  were  performed  with  linear  basis).  The  integrals  (7)  are 
then  evaluated  using  Gauss  quadrature  resulting  in  a  system  of  non-linear  ordinary  differential 
equations  which  can  be  written  as 

MY,t  =  R(Y)  (9) 

where  the  check  is  added  to  make  clear  that  Y  is  the  vector  of  solution  values  at  discrete  points 
(interpolated  through  the  space  by  the  basis  functions)  and  Yit  are  the  time  derivative  values 
at  the  same  points.  An  implicit,  second  order  accurate  family  of  time  integrators  has  been 
developed  [16]  for  application  to  this  problem.  This  new  time  integrator  is  an  extension  of 
Chung  and  Hulbert’s  generalized  alpha  method  [17]  to  our  first  order  system.  The  family  is 
characterized  by  one  free  parameter,  the  amplification  factor  as  the  time  step  over  the  period 
of  the  frequency  of  interest  tends  to  infinity,  Poo,  which  may  be  chosen  to  have  a  value  between 
0  (Gears  Method)  and  1  (Midpoint  Rule).  Full  documentation  of  Chung  and  Hulbert’s  method 
extended  to  the  Navier-Stokes  equations  is  not  possible  in  the  confines  of  this  report 

but  should  appear  in  the  literature  shortly  [16].  The  method  results  in  a  non-linear  matrix 
problem  that  is  solved  in  a  predictor-corrector  format  yielding  successive  linear  problems.  Each 
linear  problem  is  solved  using  the  Matrix-Free  Generalized  Minimal  RESidual  (MF-GMRES) 
solution  technique  with  a  block  diagonal  preconditioner  developed  by  Johan  et  al.  [18].  Con¬ 
vergence  of  the  non-linear  problem  is  confirmed  before  moving  to  the  next  time  step.  We  have 
observed  that  there  is  little  difference  in  the  turbulent  statistics  between  solving  the  equations 
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with  two  iterations  (2. 5-3.0  order  of  magnitude  reduction)  and  three  iterations  (3. 0-4.0  orders 
of  magnitude  reduction).  There  is  of  course  a  difference  in  the  instantaneous  signal  after  a  long 
time  integration  as  will  always  be  the  case  when  dealing  with  turbulence.  It  is  interesting  to 
note  that  even  roundoff  errors  will  eventually  lead  to  a  divergence  of  an  instantaneous  signal 
regardless  of  the  level  of  solution.  This  is  initially  alarming  until  we  remember  that  it  is  the 
turbulence  statistics  that  we  are  most  interested  in  predicting. 


2.2  Parallel  implementations 

The  original  implementation  of  this  approach  was  on  a  Thinking  Machines  CM5  following  the 
so-called  data  parallel  approach.  Finite  element  methods  are  amenable  to  this  approach  since 
the  bulk  of  the  computational  effort  lies  in  evaluating  local  integrals  over  each  element  domain. 
Parallelization  of  these  operations  is  trivial  after  a  gather  of  the  data  from  the  global  nodes  to 
an  element  based  data  structure.  After  the  element  integrations  are  performed  in  parallel,  the 
results  (residuals  of  (9))  are  then  scattered  (assembled)  back  to  the  nodes  where  the  equations 
are  solved. 

On  the  CM5  these  gather /scatter  operations  are  performed  using  CMSSL  library  procedures 
requiring  no  pre-processing  of  the  data.  This  made  application  to  this  platform  quite  easy. 
Unfortunately,  other  parallel  platforms  do  not  have  libraries  of  this  type.  For  this  reason  a 
second  version  of  the  code  was  developed,  employing  the  Message  Passing  Interface  (MPI) 
library  to  do  the  gather  and  scatter  operations.  This  version  of  the  code  requires  that  the 
problem  be  pre-processed  to  break  the  computational  domain  into  n  pieces  where  n  is  the 
number  of  processors,  on  which,  the  problem  will  be  run.  This  pre-processing  step  is  non-trivial 
for  unstructured  grids  since  load  balance  (equal  distribution  of  the  elements  onto  each  processor) 
and  minimal  communication  are  necessary.  Recently  we  have  accomplished  this  pre-processing 
procedure  and  have  begun  testing  our  MPI  version  of  the  code  on  a  variety  of  platforms  (IBM- 
SP2,  SGI-Origin,  SUN-Sparc  and  Cray  J90).  The  algorithm  has  shown  perfect  scalability  (see 
Bastin  [19])  on  moderate  sized  problems  (typically  the  bigger  the  better)  due  to  the  low  ratio 
of  communication  to  calculation  that  is  inherent  in  finite  element  methods.  The  MPI  version 
of  the  code  also  performs  the  communication  much  more  efficiently  than  the  CMSSL  version  of 
the  code.  On  the  CM5,  communication  required  20-25%  of  the  total  CPU  time.  With  the  MPI 
code  this  percentage  is  reduced  to  3-5%.  While  this  improvement  is  dramatic,  both  methods 
still  use  the  majority  of  the  CPU  time  computing,  making  communication  improvement  have 
little  impact.  Greater  efficiency  gains  require  a  reduction  in  the  computational  effort,  even  if  it 
comes  at  some  loss  in  communication  efficiency. 

2.3  Dynamic  model  development 

The  full  development  of  the  dynamic  model  requires  more  space  than  is  available  here.  We  refer 
the  interested  reader  to  [9],  Moin  &  Jimenez  [10],  Germano  et  al.  [11].  Here  we  only  summarize 
our  implementation  and  improvements 

Most  LES  models  rest  on  the  assumption  that  the  subgrid-scale  stresses  can  be  modeled  by 

4  =  -2/*r  Sfj(u)  (10) 
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where  the  eddy  viscosity,  fix ,  is  usually  given  by  the  Smagorinsky  model  [20], 

ttr  =  CA2?|S|,  |S|  =  fajSii  (11) 

where  C  is  the  so-called  Smagorinsky  “constant” . 

This  “constant”  proved  very  elusive  to  the  researchers  of  the  seventies  and  the  eighties. 
After  several  studies  it  seemed  clear  that  this  “constant”  in  fact  varied  from  flow  to  flow  and, 
worse  yet,  seemed  to  be  a  function  of  position  within  even  the  simplest  flows  (especially  in 
near  wall  regions).  In  the  early  nineties  a  dynamic  procedure  for  determining  the  subgrid-scale 
viscosity  was  developed  by  Germano  et  al.  [11].  In  this  procedure  the  value  of  the  constant 
C  is  determined  dynamically  by  comparing  the  solution  and  various  functions  of  the  solution 
obtained  from  the  simulation  to  the  same  quantities  passed  through  a  filter 

Applications  of  the  Smagorinsky-based  dynamic  model  flourished  within  the  spectral  meth¬ 
ods  community  and  the  finite  difference  community  where  it  gained  wide  acceptance.  Use  of 
this  model  on  unstructured  grids  required  an  extension  of  the  filtering  operator.  A  variety  of 
filtering  operators  were  developed  in  [1].  Through  our  recent  research  ([5],  [3],  [2]),  we  have 
determined  that  the  generalized  top-hat  filter  represents  the  best  balance  between  accuracy  and 
efficiency. 

2.4  Numerical  simulations  and  validations 

Reference  publications  [16,  21,  22,  23,  2,  24,  3,  25,  1,  26,  4,  5,  27]  provide  a  full  description  of 
our  work  in  this  area  but  here  we  will  summarize  it  as: 

1.  With  the  dynamic  model  disabled  we  used  our  code  to  predict  the  growth  of  the  most 
unstable  eigenmode  of  the  Orr-Sommerfield  equations  [3].  In  that  reference  we  showed 
that  our  method  is  significantly  more  accurate  than  staggered-grid  central  differences,  the 
workhorse  of  the  LES  community.  Then,  with  the  model  turned  back  on  we  confirmed 
that  the  value  of  C  vanished  when  the  disturbance  was  adequately  represented. 

2.  The  dynamic  model  was  validated  by  studying  the  decay  of  isotropic  turbulence,  where 
good  prediction  both  of  the  decay  rate  and  the  shape  of  the  spectra  were  observed  [3]. 

3.  The  classic  case  of  turbulent  flow  in  a  channel  at  Rer  =  180  was  studied  and  compared  to 
staggered-grid  central  difference  results  [3].  Again,  superior  accuracy  was  demonstrated. 

4.  Moving  to  more  complex  geometry,  a  turbulent  flow  over  a  cavity  is  currently  being 
simulated  (Reg  =  3500,  Ren  =  27000,  M  =  0.1).  Results  have  not  been  published  but  an 
instantaneous  visualization  of  the  streamwise  vorticity  is  included  here. 

5.  The  airfoil  simulation  has  been  computed  and  is  described  in  several  references  ([1],[4], 
[5], [3]).  An  instantaneous  visualization  of  the  isosurfaces  of  the  spanwise  velocity  is  in¬ 
cluded  in  Figure  2.  The  most  challenging  aspect  of  this  flow  is  the  prediction  of  smooth 
separation  on  the  last  20%  of  the  upper  surface.  The  early  LES  predictions  showed  a 
larger  separation  than  the  RANSS  simulations  but  still  less  than  that  of  the  experiments 
([28], [29], [30]).  Some  insight  into  the  reason  for  the  disagreement  can  be  obtained  from 
Figure  3.  Here,  the  calculation  is  in  good  agreement  with  the  experimental  data  when 
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Turbulent  flow  over  a  cavity 


Figure  1:  Isosurfaces  of  vorticity  in  a  turbulent  boundary  layer  flowing  over  a  cavity. 
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Figure  2:  Tip-to-tail  development  of  the  turbulent  structures  identified  by  spanwise  velocity  isosurfaces  in  the  boundary  layer 
over  the  upper  surface  of  the  airfoil.  Top  figure  provides  a  side  view  which  illustrates  the  development  of  the  boundary  layer 
thickness.  The  middle  figure  tilts  slightly  to  show  the  spanwise  structures.  Lower  left  figure  zooms  in  on  the  early  boundary 
layer  while  the  lower  right  zooms  on  the  tail  region  illustrating  the  dramatic  variation  in  spatial  scales. 


Figure  3:  Coefficient  of  pressure  on  the  airfoil  surface  comparing  calculation  ( - )  to  Hasting’s 

experiment  with  transition  strip  (o  )  and  without  transition  strip  (+  ). 


there  is  no  transition  strip.  However,  there  is  a  substantial  flattening  of  the  experimental 
Cp  distribution  in  the  trailing  edge  region  when  the  transition  strip  was  in  place.  This 
minor  change  in  the  Cp  distribution  signals  a  much  more  dramatic  change  in  the  velocity 
field  where  a  much  more  massive  separation  would  be  expected.  While  it  is  encouraging 
that  our  calculation  (which  at  this  point  had  no  transition  strip)  agrees  with  the  experi¬ 
ment  with  the  same  boundary  conditions,  all  velocity  and  stress  data  was  taken  with  the 
transition  strip  in  place.  Therefore,  in  subsequent  calculations  we  modified  our  boundary 
conditions  to  accurately  represent  both  the  wind  tunnel  walls  and  the  serrated  transi¬ 
tion  strip  of  Wadcock  [30]  as  seen  if  Figure  4.  Addition  of  these  features  did  improve 
the  agreement,  but  discrepancies  remained.  We  believe  that  the  remaining  discrepancies 
are  linked  to  the  inadequate  spanwise  extent  of  our  computational  domain.  Figure  5 
shows  the  effect  of  doubling  the  spanwise  domain  width  on  the  velocity  profiles  in  the 
separated  region.  The  obvious  solution  of  further  doubling  the  spanwise  domain  width  is 
unattractive  since  it  would  again  double  the  cost  of  the  calculation.  Instead,  under  other 
support,  we  are  developing  new  boundary  conditions  to  address  this  difficulty  in  a  more 
cost  effective  manner. 

6.  The  above  computations  were  completed  using  only  linear  interpolation  functions.  We 
have  recently  implemented  a  new  set  of  interpolation  functions  which  utilize  a  hierarchical 
basis  [21,  22].  With  one  input  parameter  change,  we  can  change  the  order  of  accuracy  of 
our  method.  We  have  confirmed  that  we  obtain  the  theoretically  optimal  rate  of  conver¬ 
gence  ( 0(hp+1 )  in  the  L2  norm  where  p  is  the  polynomial  order)  for  laminar  channel  flow 
and  acoustic  wave  propagation  and  for  Kovasnay  flow.  This  approach  adds  p  adaptivity 
to  our  already  demonstrated  capability  in  h  adaptivity.  While  other  researchers  are  also 
exploring  this  approach  we  should  emphasize  that  we  have  successfully  integrated  this 
approach  into  our  massively  parallel,  dynamic  subgrid-scale  model  LES  code.  Therefore, 
we  have  removed  all  of  the  obstacles  to  large  scale  simulations  with  this  approach.  Indeed, 
turbulent  simulations  are  underway  at  the  time  of  this  writing  and  preliminary  results 
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Figure  4:  Geometric  model  of  the  exact  shape,  height  and  position  (shown  in  black  on  the 
airfoil  surface)  of  Wadcock’s  transition  strip. 


Figure  5:  Profiles  of  tangential  velocity  component  at  various  positions  along  the  airfoil  surface 
(x/c  =  0.59,0.066,0.78,0.82,0.95,  plots  have  been  shifted  by  1.5  at  each  station).  Solutions 

correspond  to:  without  wind  tunnel  walls  or  transition  strip - ,  with  wind  tunnel  walls  and 

transition  strip  ( W/c  —  0.025) - ,  with  wind  tunnel  walls  and  transition  strip  (W/c  =  0.05) 

- ,  Wadcock  o  ,  Hastings  and  Williams  +  . 


are  quite  promising.  We  expect  to  present  them  at  the  Second  AFOSR  International 
Conference  on  DNS/LES. 

2.4.1  Improved  solver  efficiency 

We  currently  use  a  matrix-free  generalized  minimal  residual  method  (MF-GMRES),  which 
requires  very  low  storage  but  is  somewhat  computationally  intensive.  Since  memory  is  becoming 
more  available  on  parallel  machines  we  have  investigated  a  new  solver  which  takes  advantage 
of  the  sparseness  of  our  matrix,  thereby  enabling  us  to  form  the  matrix  which  greatly  reduces 
the  computational  effort  at  a  minor  cost  to  memory.  Our  current  experience  with  the  8.5 
million  element  airfoil  simulation  [2]  suggests  that  this  tradeoff  will  be  a  highly  favorable  one. 
Forming  the  matrix  has  a  second  benefit  since  it  allows  a  much  richer  set  of  pre-conditioners 
to  be  employed  (currently  in  the  matrix-free  technique,  only  block-diagonal  pre-conditioning  is 
possible),  further  accelerating  convergence  at  each  step.  AFOSR  support  produced  a  Masters 
thesis  on  this  topic  which  is  in  the  process  of  being  converted  into  a  journal  article  (thesis  is 
available  by  request). 

2.4.2  New  boundary  conditions 

A  very  promising  avenue  to  improved  simulations  is  through  the  introduction  of  new  boundary 
conditions.  Here  we  added  two  new  types  of  boundary  conditions  to  our  method.  The  first, 
a  sub-layer  model,  is  somewhat  classical.  The  second,  and  perhaps  more  novel,  can  best  be 
described  as  a  surface  extraction  boundary  condition  (SEBC).  In  the  following  paragraphs  we 
describe  each  of  these. 

Wall  layer  modeling  has  been  carried  out  since  the  early  days  of  LES  [31].  It  is  known  to  be  a 
fairly  good  approximation  for  flows  that  are  near  equilibrium  but  can  give  disappointing  results 
in  other  cases.  The  advent  of  unstructured-grid  LES  affords  the  opportunity  to  develop  new 
wall  layer  models  which,  through  their  dependence  on  grid  density  could  smoothly  fair  between 
a  resolved- wall  LES  and  a  wall-layer-modeled  LES.  In  this  way,  a  full  LES  could  be  carried  out 
in  regions  where  it  is  necessary  (i.e.  where  the  flow  is  not  in  equilibrium  and  high  accuracy 
is  needed)  and  a  wall  layer  model  could  be  used  elsewhere  (where  the  flow  is  in  equilibrium 
and/or  accuracy  requirements  are  lower).  Existing  wall  layer  models  do  not  have  this  character 
and  therefore  new  models  must  be  developed.  Again,  we  emphasize  that  the  unstructured-grid 
LES  approach  offers:  1)  the  opportunity  to  resolve  more  interesting  cases  in  the  validation  of 
these  models,  and  2)  the  capability  to  smoothly  change  the  mesh  spacing  between  the  fine  near 
wall  spacing  required  in  a  resolved-wall  LES  and  the  coarse  near  wall  resolution  required  for  a 
wall-layer-modeled  LES.  We  have  completed  preliminary  testing  of  the  wall  layer  model  wherein 
we  have  applied  it  to  channel  and  flat  plate  flows.  Application  to  the  airfoil  is  the  next  step. 

Experience  with  the  airfoil  LES  has  identified  a  new  boundary  condition  that  will  dramat¬ 
ically  reduce  the  computational  effort  required  in  LES  of  problems  with  a  solid  surface  whose 
dimension  transverse  to  the  mean  flow  direction  is  large  (e.g.  the  leading  edge  of  a  wing).  In 
problems  of  this  type,  what  makes  LES  so  expensive  is  the  need  to  carry  the  fine  resolution  to 
resolve  the  turbulent  fluctuations  over  a  long  transverse  or  spanwise  dimension.  To  accurately 
model  the  turbulent  fluctuations  in  this  region  only  requires  a  very  short  spanwise  dimension 
since  they  are  known  to  have  a  fairly  short  two  point  correlation  decay  length  (typically  600-800 
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wall  units  of  spanwise  extent  are  necessary  to  insure  a  de-correlated  signal).  That  is  to  say,  in 
the  leading  edge  region,  the  same  solution  statistics  would  be  obtained  for  a  spanwise  domain  of 
700  wall  units,  7000  wall  units,  or  even  70,000  wall  units.  Indeed  this  is  the  very  principle  upon 
which  the  spanwise  periodic  boundary  condition  is  based.  Typically  we  choose  our  spanwise 
domain  length  to  be  greater  than  the  two  point  correlation  decay  length  at  all  points. 

In  the  airfoil  simulation  the  two  point  correlation  length  changes  by  over  an  order  of  mag¬ 
nitude  due  to  the  strong  adverse  pressure  gradient  which  drives  the  flow  to  separate  in  the 
last  20%  of  the  chord  length.  Therefore,  our  current  simulation  is  about  8  times  wider  than  is 
necessary  for  the  early  turbulent  boundary  layer  and  perhaps  is  still  too  narrow  in  the  trailing 
edge  region.  Indeed,  when  we  compare  the  statistics  of  this  solution  to  the  statistics  of  the 
solution  obtained  on  a  grid  with  the  exact  same  grid  spacing  but  with  half  the  spanwise  length 
(translating  into  a  factor  of  two  savings  in  grid  points),  the  early  boundary  layer  is  unchanged 
while  the  trailing  edge  boundary  layer  is  adversely  impacted.  By  examining  Figure  2  one  can 
see  that  the  solution  near  the  nose  looks  to  be  de-correlated  in  about  one  eighth  of  the  domain 
(i.e.  the  typical  eddy  structure  is  about  one  sixteenth  of  the  spanwise  length),  while  the  trailing 
edge  boundary  layer  contains  only  one  or  two  pairs  of  structures. 

We  have  completed  preliminary  testing  of  a  new  boundary  condition  to  alleviate  this  diffi¬ 
culty.  This  boundary  condition  allows  the  domain  width  to  be  doubled  as  necessary  to  maintain 
an  adequate  two  point  correlation  decay  length  while  maintaining  computational  efficiency.  To 
correctly  implement  this  boundary  condition  requires  that  we  find  the  appropriate  solution  to 
apply  as  a  boundary  condition  on  the  newly  doubled  surface.  The  appropriate  solution  is  the 
solution  from  the  plane  that  was  just  copied.  The  procedure  is  illustrated  for  the  surface  nodes 
in  Figure  6.  The  process  must  be  carried  out  for  the  plane  extending  from  these  nodes  to  the 
opposing  boundary  but  for  clarity  we  only  show  the  line  here.  We  name  this  boundary  condition 
a  surface  extraction  boundary  condition  (SEBC).  The  SEBC  is  effectively  an  inflow  boundary 
condition  where  the  solution  value  is  set  to  that  of  a  companion  node  on  the  interior  plane. 
The  process  can  be  repeated  as  many  times  as  is  necessary  with  the  only  limitation  being  that 
the  spanwise  width  must  grow  by  an  integer  multiple  of  the  previous  spanwise  width. 

Though  conceptually  simple,  there  remain  a  number  of  challenges  to  realizing  the  efficiency 
gains  of  the  SEBC.  The  key  points  that  must  be  addressed  in  this  research  are: 

1.  An  initial  mesh  is  generated  that  is  wide  enough  for  the  most  narrow  section, 

2.  The  appropriate  points  to  expand  the  domain  are  identified  and  the  mesh  is  copied  to 
create  a  discretely  varying  spanwise  width  (all  of  the  nodes  and  elements  between  the  red 
plane  and  the  opposite  boundary  in  Figure  6), 

3.  At  the  same  time,  the  newly  created  nodes  which  fall  on  the  inflow  boundary  are  given  a 
boundary  condition  which  links  them  to  the  nodes  that  they  were  copied  from  (the  nodes 
along  the  plane  created  by  extending  the  purple  line  to  the  opposite  boundary  are  linked 
to  the  nodes  between  the  green  line  and  the  opposite  boundary), 

4.  The  other  newly  created  nodes  which  fall  on  the  other  boundaries  are  given  the  appropriate 
boundary  conditions,  (including  the  periodic  boundary  condition  which  is  transferred  from 
the  points  which  are  now  interior  to  the  newly  expanded  cross  sectional  plane), 
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domain  is  doubled  (shown  in  red).  The  inflow  to  the  new  domain  (the  plane  between  the  purple  line  and  the  farfield  boundary) 
is  linked  a  plane  on  the  interior  of  the  original  domain  (the  plane  between  the  green  line  and  the  farfield  boundary).  The  process 
is  repeated  again  at  a  later  section  through  a  second  doubling  of  the  spanwise  domain.  We  only  show  the  surfaces  here  for 
clarity  but  we  emphasize  that  the  boundary  condition  applies  to  the  newly  created  inflow  planes. 


5.  Though  this  boundary  condition  does  not  initialize  any  energy  into  the  new  low  wavenum¬ 
ber  modes  that  are  supported  on  the  expanded  domain,  this  is  a  spatially  developing  flow 
and  these  modes  can  be  expected  to  grow  in  the  new  domain  where  they  are  no  longer 
constrained  by  the  periodic  boundary  condition.  Preliminary  testing  on  flate  plate  and 
channel  flows  has  been  completed  with  encouraging  results.  Application  to  more  complex 
geomtries  is  underway. 

As  noted  in  [32],  it  is  the  boundary  layer  near  the  leading  edge  which  has  the  most  strict 
resolution  requirements  leading  to  the  prohibitively  large  number  of  grid  points  required  for 
wing  geometries.  This  approach  is  expected  to  reduce  the  current  airfoil  simulation  cost  by  an 
order  of  magnitude.  The  wing  simulation  cost  is  expected  to  drop  by  two  to  three  orders  of 
magnitude.  These  gains  are  possible  because  the  leading  edge  region  is  the  most  likely  region 
to  satisfy  the  conditions  mentioned  above  for  application  to  spanwise  inhomogeneous  flows. 

2.5  Project  summary 

The  work  summarized  herein  aggressively  pushed  the  envelope  of  turbulence  simulation  through 
the  development  of  more  efficient  methods,  solvers,  models,  and  boundary  conditions. 
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